
#include "VAR.hpp"

//
// constructor
// initialize everything
//
VAR::VAR ( int p_order, int p_dim ) :
  m_order(p_order), m_dim(p_dim), m_A(p_order), m_Y(p_order,p_dim,0.0),
  m_Yt(p_dim,0.0), m_C(p_dim,0.0), m_S(p_dim,1.0), m_E(p_dim,0.0)
{
  RN = new Normal;
  // the coefficient matrices have not been initialized yet
}

VAR::VAR ( int p_order, int p_dim, const vector< Matrix<double> >& p_A,
	   const Vector<double>& p_C ) :
  m_order(p_order), m_dim(p_dim), m_A(p_A), m_C(p_C),
  m_Y(p_order,p_dim,0.0), m_Yt(p_dim,0.0), 
  m_S(p_dim,1.0), m_E(p_dim,0.0)
{
  RN = new Normal;

  if ( p_A.size() != m_order )
    throw (range_error("VAR::VAR: matrix size (order) mismatch"));
  if ( (p_A[0]).size() != m_dim )
    throw (range_error("VAR::VAR: matrix size (dimension) mismatch"));
}

VAR::~VAR()
{
  delete RN;
}

void
VAR::Sigma ( const Vector<double>& p_S )
{
  if ( p_S.size() != m_dim )
    throw (range_error("VAR::Sigma: size mismatch."));

  m_S = p_S;
}

void
VAR::Yi ( const Vector<double>& p_Yi, int i )
{
  if ( i >= m_order )
    throw (range_error("VAR::Yi(Vector,int): subscript out of range."));
  m_Y[i] = p_Yi;
}
const Vector<double>&
VAR::Yi ( int i )
{
  if ( i >= m_order )
    throw (range_error("VAR::Yi(int): subscript out of range."));
  return m_Y[i];
}

void
VAR::Ai ( const Matrix<double>& p_Ai, int i )
{
  if ( i >= m_order )
    throw (range_error("VAR::Ai(Matrix,int): subscript out of range."));
  m_A[i] = p_Ai;
}
const Matrix<double>&
VAR::Ai ( int i ) const
{
  if ( i >= m_order )
    throw (range_error("VAR::Ai(int): subscript out of range."));
  return m_A[i];
}
void
VAR::C ( const Vector<double>& p_C )
{
  m_C = p_C;
}

//
// one step forward
//
const Vector<double>&
VAR::simulate ()
{
  // generate a new disturbance
  for ( int i=0; i < m_dim; ++i )
    m_E[i] = RN->rand() * m_S[i];

  // add constant component
  m_Yt = m_C + m_E;

  for ( int i=0; i < m_order; ++i )
    m_Yt += m_A[i] * m_Y[i];

  // shuffle data
  // the last piece of data is lost
  for ( int i = m_order-1; i > 0; --i )
    m_Y[i] = m_Y[i-1];
  m_Y[0] = m_Yt;

  return m_Yt;
}

